Laser harmonic generation with independent control of frequency and orbital angular momentum

The non-linear optical process of laser harmonic generation (HG) enables the creation of high quality pulses of UV or even X-ray radiation, which have many potential uses at the frontiers of experimental science, ranging from lensless microscopy to ultrafast metrology and chiral science. Although many of the promising applications are enabled by generating harmonic modes with orbital angular momentum (OAM), independent control of the harmonic frequency and OAM level remains elusive. Here we show, through a theoretical approach, validated with 3D simulations, how unique 2-D harmonic progressions can be obtained, with both frequency and OAM level tuned independently, from tailored structured targets in both reflective and transmissive configurations. Through preferential selection of a subset of harmonic modes with a specific OAM value, a controlled frequency comb of circularly polarised harmonics can be produced. Our approach to describe HG, which simplifies both the theoretical predictions and the analysis of the harmonic spectrum, is directly applicable across the full range of HG mechanisms and can be readily applied to investigations of OAM harmonics in other processes, such as OAM cascades in Raman amplification, or the analysis of harmonic progressions in nonlinear optics.

In this section, we expand upon the technique used throughout the main text to diagnose the harmonic spectrum obtained from our 3D simulations.We are obliged to study the electric field ⃗ E, since our simulation results do not include the vector potential ⃗ A. A transverse real-valued electric field ⃗ E ⊥ can be decomposed into two orthogonal linear polarisations: ⃗ E ⊥ = ( ⃗ E ⊥ •⃗ e x )⃗ e x + ( ⃗ E ⊥ •⃗ e y )⃗ e y , where ⃗ E ⊥ •⃗ e x and ⃗ E ⊥ •⃗ e y are independent real-valued quantities.Alternatively, ⃗ E ⊥ can also be decomposed into two orthogonal circular polarisations: ⃗ E ⊥ = E + ⃗ e + + E − ⃗ e − , with ⃗ e ± = (⃗ e x ± i⃗ e y )/ √ 2 and E ± = ⃗ E ⊥ • ⃗ e ± complex-valued.However, since ⃗ E * ⊥ = ⃗ E ⊥ while ⃗ e * − = ⃗ e + , we find that E − = E * + .In the frequency domain, we find that the Fourier transform F[E − ](ω) = (F[E + ](−ω)) * , i.e. the RCP frequency spectrum for ω > 0 is the complex conjugate of the LCP frequency spectrum for ω < 0 and vice versa.This motivates us to only use the E + coefficient, or actually Ψ ≡ E x + iE y = E + √ 2 = −∂Ψ/∂t, and to use the "signed frequency" −∞ < ω/σ < ∞, as in Section S2.Thus, F[ Ψ](ω/σ) captures all the spectral information of ⃗ E ⊥ , while each value of ω/σ represents a mode with pure circular polarisation.Naturally, the use of ω/σ in the Fourier domain implies the use of ⃗ k/σ, ℓ/σ, etc. also.
This technique is especially important in the case of a circularly polarised pump and a harmonic with opposite spin to the pump, since such harmonics tend to be drowned out by their counterpart with the same spin as the pump.See, for example, Li et al. [3], in which harmonics with opposite spin are predicted but not observed, as their amplitude is assumed to be "too small" to be detected easily.

S2. CONVENTIONS
In the main text, we study laser harmonic generation (HG) by a cubic nonlinear term | ⃗ A| 2 ⃗ A in the wave equa-tion for the vector potential ⃗ A of the electromagnetic (laser) field.This covers the most common scenarios for harmonic generation.Throughout, we use the following conventions: 1. We write all laser fields as superpositions of circularly polarised, CP, laser modes.For example, a laser pulse with linear (elliptic) polarisation is the superposition of two CP modes with equal (unequal) amplitudes and opposite spin.
For two such modes ⃗ A A and ⃗ A B , we find that The nonlinear term in the wave equation for Ψ then becomes (Ψ * Ψ)Ψ.This convention will also allow us to make use of the complex Mathieu equation [4].
3. In light of the above, we use the quotients ω/σ, ⃗ k/σ, ℓ/σ, etc. instead of σ and ω separately.While ω/σ can be any real number, we always use ω > 0, which uniquely defines σ.This convention renders both the calculation and the visualisation of harmonic progressions easier and more intuitive.We use ω/σ rather than σω since ω/σ = ℏω/(ℏσ) = E/S z , with E and S z denoting the photon energy and spin.The quantity σω does not have such an elegant physical interpretation.
4. We treat left and right CP modes at the same harmonic frequency ω as separate harmonics with opposite signs of ω/σ rather than different polarisations of the same harmonic.We do this so that we can treat the spin σ on a footing equal to the frequency ω, wave vector ⃗ k or OAM level ℓ.For example, sending a left CP laser beam through a halfwave plate to create a right CP beam will be treated as generating harmonic −1, since ω/σ → −ω/σ in that case.
5. In terms of the harmonic progression, we treat harmonic cascades in quantities like ω, k ⊥ or ℓ on an equal footing as much as possible.So a cascade in ℓ at a single frequency ω 0 [5] is treated as a harmonic progression just like a cascade in ω at the single OAM level ℓ = 0.This way, multi-dimensional cascades in, for example, (ω/σ, ℓ/σ) space are also possible.
The justification for these conventions is as follows: (i) A CP wave has constant amplitude A 2 and does not beat with itself (unlike an LP wave), while two CP waves with parallel ⃗ k beat in precisely one way.
(ii) A CP wave is a spin eigenstate, while LP and EP waves are not.Thus, CP waves can be used to verify the spin conservation law, while LP waves cannot be used for that.
(iii) Since a CP wave has constant A 2 , it obeys all possible translation/rotation symmetries.Adding a second CP wave will break the symmetry of A 2 in one specific way, leading to harmonic generation in one specific dimension in (ω/σ, ⃗ k/σ) space.This relates the dimensionality of symmetry breaking to the dimensionality of the harmonic progression.
Thus, CP waves are more fundamental "building blocks" to study HG than LP waves.Using the above conventions renders it easier to calculate and visualise the harmonic progression.The harmonic progression becomes a "beat wave" process [6,7]: the beating fixes the harmonic step to ω B /σ B − ω A /σ A or ℓ B /σ B − ℓ A /σ A etc., while ω/σ (or ℓ/σ etc.) can assume any positive or negative value.It is easy to extend this procedure to more variables (k ⊥ , ℓ, etc.) or to more CP modes.The harmonic spectrum generated by two CP modes will be a sequence of equidistant points on a line; for three CP modes that are not on the same line in 2-D Fourier space, the generated harmonic spectrum will be a 2-D grid; for four CP modes that are not on the same plane in 3-D Fourier space, the generated harmonic spectrum will be a 3-D lattice.

S3. HARMONIC PROGRESSION
The details of the harmonic progression for HG in low-order nonlinear systems are outlined in this section.In most practical cases, the lowest-order nonlinear term in the wave propagation equation is either cubic, e.g.(Ψ * Ψ)Ψ, or quadratic, e.g.2(Ψ * DC Ψ)Ψ or similar, where Ψ DC is constant in time.In general, the cubic case is seen in HG in isotropic media where no energy or momentum is left behind in the medium, so energy and momentum are conserved in the EM waves alone.The quadratic case is often seen in HG in anisotropic media where linear or angular momentum (or even energy, as in stimulated Raman scattering) can be left behind in the medium.In the quadratic case, the nonlinear terms can often be rewritten as effectively reducing it to a cubic case providing a joint description of the laser and the medium, see also Ref. [8].For this reason, we first study the case of a wave equation with cubic nonlinearity: etc.For a single mode Ψ 0 = Ψ A , we find that Ψ * 0 Ψ 0 is constant, and no harmonics will be generated.For two modes, Ψ 0 = Ψ A + Ψ B , we can use the theory of the Mathieu equation and the Jacobi-Anger expansion to find the harmonic frequency progression: Since ω n > 0 always, σ n is well-defined.We note that all the frequencies in the ω/σ spectrum are equidistant, while this need not be true at all for the ω spectrum.Now that σ n has been determined, we introduce the other harmonic progressions (n ∈ Z): where the last equation is of course trivially true, and ensures that every harmonic step is spin-neutral.We define # A , # B , etc. as the number of photons absorbed from mode A, B, etc. for the generation of mode n, and make use of the fact that any "fundamental mode" is an eigenfunction for energy with eigenvalue ℏω A etc. Since the nonlinear factor Ψ * Ψ in (E1) obeys the same symmetries as the wave energy ∝ E 2 , it is safe to assume energy conservation [9][10][11][12], i.e. ω n = # A ω A + # B ω B .Then we find: We note that every quantity that scales in the same way as ω n , i.e. σ n , ℓ n , and k ⊥,n , will be automatically conserved when ω n is conserved, due to the form of the expressions for # A and # B , and due to the fact that Ψ A is an eigenfunction for all these quantities.Due to our conventions, we can calculate # A and # B explicitly, and get conservation laws for spin, orbital angular momentum and transverse momentum at the same time.
The harmonic progression (E2)-( E6) is characterised by one index n and is essentially one-dimensional.However, with three independent CP modes, a twodimensional progression in two variables, e.g.ω and ℓ, can be realized.The harmonic progression then becomes: The most visible difference between the 1-D and 2-D progressions is that in a 1-D progression there will always be exactly one OAM level for a given harmonic frequency, while in a 2-D progression there can often be many OAM levels for a given harmonic frequency.

S4. MOVING WINDOW
In order to observe the 2D harmonic spectrum after the laser pulse has propagated into the far-field after interaction with the aperture targets, additional simulations using a moving window were conducted.These were defined with a larger simulation domain of −20 µm ≤ x, y ≤ 20 µm and −25 µm ≤ z ≤ −5 µm with 1440×1440×1500 grid cells.The target and laser parameters were the same as the previous simulations with the exception that the temporal profile was reduced to 30 fs FWHM to ensure the entire pulse fit in the simulation domain.The simulation allowed the laser pulse to propagate for 60 fs before the simulation domain begins to move in z -direction with a velocity of c.This allows the entire pulse to be observed as it moves through the aperture and beyond without having to continually simulate the plasma.The pulse was then sampled when the simulation domain was centred at z =250 µm.Due to the nature of the simulation and as the pulse is now in free-space, a spatial Fourier transform is used instead of temporal and the 2D harmonic spectrum is constructed.This can be seen in Supplementary figure S1a and b for n=0 and n=3, respectively.As can be observed, the harmonic spectra maintain the same profile, albeit with more numerical noise potentially due to numerical dispersion when considering the pulse in the spatial domain instead of frequency.To ensure we are indeed in the far-field, the spatial profile of the -ω and +ω components of the generated 2ω 0 harmonic light are shown in Supplementary figure S1c and d for the n=0 case and Supplementary figure S1e and f for the n=3 case, respectively.As can be seen, these correspond to the constructed k-space patterns observed in figure 3e-h of the main text.

S5. TARGET STRUCTURE SIZE EFFECTS
The variation of the radius and corrugation depth of the targets was investigated for the transmissive cases.Supplementary figures S2 a and b show the generated 2D harmonic spectra for an aperture radius of 1.5 µm and corrugation depths of 0.5 µm and 1.5 µm, respectively.The overall structure of the spectra is relatively unchanged, albeit the strengths of the different harmonic peaks vary.The same can be seen in Supplementary figures S2c-d for the same corrugation depths, but an increased radius of 2.5 µm.The change in strength of the harmonics can impact on the overall spatial profile of the generated harmonics as shown in the full 2ω 0 2D profile shown Supplementary figures S2e-h for the same preceding cases, respectively.

S6. FROM PLASMA PHYSICS TO NONLINEAR OPTICS
Our approach provides access to a vast landscape of extensions, both in laser-plasma interactions and in nonlinear optics.Within plasma physics we find: (i) Adding OAM to the pump laser, changes the shape of the 2-D grid and thus the frequencies in the comb.
(ii) A tilted target adds a DC mode with ℓ/σ = 0, while more complex targets can add multiple values of ℓ/σ.
(iii) Qu, Jia and Fisch [18] studied a plasma q-plate for q = 1 using the magnetic field of two Helmholtz coils.We can now study plasma q-plates for q ≤ 0 using the magnetic field of a magnetic 2(1 − q) multipole (as used for particle beam optics), which already has the required shape.
(iv) By adding OAM to the pump laser, the concept of the non-linear q-plate [19,20] can be extended to that of the non-linear J-plate [21,22], provided we use ℓ/σ of both the pump laser and the target as two independent degrees of freedom.
Since HG in nonlinear optics is usually governed by the same cubic nonlinearities, many of our findings can be applied to this field without much modification.Examples include: (i) The work by Hickstein et al. [23], which reports on HG by CP beams (two fundamental modes, 1-D expansion, one transverse angle per frequency) versus LP beams (four modes, 2-D expansion, multiple transverse angles per frequency).The findings of that work can be explained easily using our approach.
(iii) The "torus-knot angular momentum", introduced in Pisanty et al. [25] and Dorney et al. [26]: the two fundamental modes generate a sequence of harmonic modes that are all on the same line in (ω/σ, ℓ/σ) space, at regular intervals, without the need to introduce new physical quantities.
(iv) In stimulated Raman scattering [27,28] or beat wave generation [6,7], a laser beam interacts with a moving structure in a nonlinear medium to produce side bands.This can be studied using the same approach as HG, including extension from frequency to OAM side bands [5].
These examples illustrate how general our new approach is, and the breadth of its field of applications.

S7. HEIGHT FUNCTION
In this section, we move beyond the Lichters model [13] and study HG by a laser pulse incident on an arbitrary surface of an opaque target.To describe the target surface, we use a height function h(x, y) and define ⃗ a DC = ∇ ⊥ h or Ψ DC = ∂h/∂x + i∂h/∂y.For oblique incidence on a plane target this approach yields h(x, y) = x tan α and Ψ DC = tan α.For more intricate surface shapes, there are almost unlimited possibilities for Ψ DC .
Wang et al. [14] reports on the interaction of a circularly polarised laser beam, σ = +1, with a conically shaped dent in the target surface, half angle α.Li et al. [3] reports on the interaction of a tightly focused laser pulse having a concave phase front with a planar target, which is sufficiently similar to Ref. [14] for our purposes.The height function for a conical dent is h(x, y) = tan(π/2 − α) x 2 + y 2 so Ψ DC = tan(π/2 − α) exp(iφ), where φ is the azimuthal angle in polar coordinates.Then the DC mode has ω/σ = 0 and ℓ/σ = −1.If this dent is hit by a circularly polarised laser beam with spin σ 0 , frequency ω 0 and OAM level ℓ 0 , then we find from (E2) and (E3): ω n /σ n = n(ω 0 /σ 0 ) and ℓ n /σ n = n(ℓ 0 /σ 0 + 1) − 1 for all n ∈ Z.These results encompass the selection rules provided in Refs.[3,14], and also include the "negative frequencies" for n < 0, which have been predicted in Li et al. [3] but not yet demonstrated.
Since a target with circular symmetry cannot absorb net angular momentum, we note that laser-target configurations with full rotational symmetry usually obey ℓ A /σ A = −1 and thus exhibit "conservation of the total angular momentum"; see e.g.Refs.[3,14,15,29,30].Conversely, a q-plate flips the polarisation of circularly polarised light and adds orbital AM to it, absorbing 2 units of AM and emitting 2q units [19,20].The total angular momentum of the beam then changes by 2(q − 1), which is nonzero for q ̸ = 1.
Combining all of the above, we wish to obtain a laserplasma interaction that generates harmonics satisfying (n ∈ Z, q ∈ Q): This generalises all of the above.Setting q = 0 returns the results in Lichters et al. [13].Setting q = 1 returns the results in Refs.[3,14].Setting n = −1 and thus σ n = −σ yields an ordinary q-plate [19,20].We note that Eqns.(E10) and (E11) correspond to Eqns.(E2) and (E3) in Section S3 with ω A = 0 and ℓ A /σ A = −q.Thus, we need to create a DC mode with Ψ DC ∝ exp(iqφ), i.e. a shaped target surface with a height function such that ∇ ⊥ h ∝ (cos(qφ), sin(qφ)).To determine a suitable h(x, y), we first identify its contour lines, i.e. the curves (r(φ), φ) in polar coordinates for which the polar angle α of the tangent to the curve satisfies α(φ) = qφ + α 0 .This leads to the following equation: For q ̸ = 1, we can set α 0 = π/2 without loss of generality, and we find (C > 0): For q = 1 and α 0 ̸ = 0, we find: For q = 1 and α 0 = 0, we find: φ = C, r undetermined.The next step is to determine the height function h(x, y) or h(r, φ) that we need in the laser-target interaction where we generate our harmonics.In Ref. [14], a cone-shaped dent is used with height function h(r, φ) = r tan(π/2 − α), to obtain laser-plasma HG with q = 1.The contour lines for this shape satisfy: z = C, r(φ) = C, i.e Eq. (E14) with q = 1 and α = π/2.This can easily be generalised.For q ̸ = 1, we set h(r, φ) = C 1−q and invert Eq. (E13) to obtain: For q > 1, this leads to a pole at r = 0, so we may need to exclude the area r < ϵ for some small ϵ > 0, or use e.g.h(r, φ) = arctan{r 1−q cos[(q − 1)φ]} to saturate h near r = 0.This is not necessarily a problem, since laser beams with nonzero OAM usually have vanishing intensity for small r, so the precise shape of the target surface near r = 0 is not too important.As an alternative, one can set h = C q−1 when q > 1 and invert Eq. (E13) to obtain h(r, φ) = r q−1 / cos[(q − 1)φ], and saturate h near the angles where cos[(q − 1)φ] = 0.For q = 1 and α = π/2, we recover the conical shape of Ref. [14].For q = 1 and α = 0, we set h(r, φ) = λ 0 φ/(2π) for r ≥ ϵ > 0 to obtain the "spiral plate" of Ref. [3].For q = 0, we obtain h(r, φ) = r cos(φ) = x, which corresponds to a laser beam with oblique incidence onto a plane, as described in Lichters et al. [13].We see that we can design target surfaces to obtain the harmonic generation patterns dictated by Eqns.(E10) and (E11) for any q.
The salient features of the non-linear q-plate are as follows: (i) We can have any n ∈ Z, not just n = −1 or n > 0.
(ii) We can have any q ∈ Q, not just q = 0 or q = 1, or (half-)integer q.
(iii) We can have q < 0 just as well as q > 0.
(iv) Surprisingly, for non-integer q, ℓ n given by Eq. (E11) need not be an integer.However, for integer ℓ n , the OAM spectrum for harmonic ω n /σ n will show a sharp peak at ℓ n /σ n , while for non-integer ℓ n the OAM spectrum will show a much broader range of peaks at integer OAM levels.This distinction might allow one to select or suppress certain harmonics (v) We note that for J = L + S, we find that j n /σ n = ℓ n /σ n + 1 = n(ℓ/σ + q) − (q − 1) = n(ℓ/σ + 1) + (n − 1)(q − 1), so J is conserved within the EM waves for all harmonics if and only if q = 1.For q ̸ = 1, total angular momentum is not conserved in the EM waves alone; the medium will have to absorb the difference.
The concept of the non-linear q-plate can be extended to that of the non-linear J-plate [21,22].This works as follows.A basic J-plate turns a circularly polarised beam with σ = +1 into one with σ = −1, adding ℓ r to the OAM level.Similarly, it turns a circularly polarised beam with σ = −1 into σ = +1, adding ℓ s to the OAM level.To ensure that ℓ r and ℓ s are fully independent, we need two degrees of freedom in the design of the non-linear q-plate, which we choose to be q and ℓ J , where ℓ J is the amount of OAM that needs to be added to the laser pulse prior to its interaction with the q-plate.We revisit Eq. (E11), with nonzero ℓ J : We set n = −1 and use σ −1 = −σ to obtain: For σ = +1, we find ℓ −1 − ℓ = ℓ r = ℓ J + 2q, while for σ = −1, we find ℓ −1 −ℓ = ℓ s = ℓ J −2q.Thus, a non-linear J-plate with parameters ℓ r , ℓ s can be produced from an enhanced non-linear q-plate by setting q = (ℓ r − ℓ s )/4 and ℓ J = (ℓ r + ℓ s )/2.Note that ∇ ⊥ h of a non-linear q-plate corresponds to the "frozen" field of a Laguerre-Gaussian mode with circular polarisation and OAM level ℓ/σ = −q, see e.g.Courtial et al. [31].We can exploit this to design target surfaces that correspond to e.g. the frozen field of Hermite-Gaussian modes.Irradiate this with a Hermite-Gaussian mode with linear polarisation parallel to ∇ ⊥ h, to obtain odd harmonics with the HG mode of the pump, even harmonics with the HG mode of the target.There are many more possibilities, e.g. also exploring the radial mode of the Laguerre-Gaussian harmonics, but that is beyond the scope of this work.

S8. TRANSMISSIVE AND REFLECTIVE LASER PARAMETERS
In the main text, the laser parameters for the reflective and transmissive cases are slightly different due to originally being conducted independently to verify the harmonic predictions.The simulations for the transmission case have also been run with the same parameters as the reflective case (I max = 6.6 × 10 18 Wcm −2 and 100 fs FWHM pulse duration) and the results for the n=0,2 and 5 cases are shown in Supplementary figure S3a-c.As can be seen, the overall harmonic spectra maintains the same profile, but the overall strength of the harmonics is reduced.This is particularly noticeable for the higher harmonics.The on-axis y = 0 and x = 0 spectra are also shown in red in Supplementary figure S3d-f with the original I max = 1 × 10 19 Wcm −2 and 50 fs FWHM pulse duration is shown in gray for comparison.For the longer pulse, lower intensity, the 6ω 0 can no longer be observed above the numerical noise.In order to maintain better signal-to-noise, the higher intensity shorter pulse duration results are used throughout in the main text for the transmission cases.

Figure
Figure S1: a-b, 2D harmonic spectrum for the, a, 0-fold and, b, 3-fold apertures from the moving window simulation.c-d, Spatial profile of the generated -2ω0 and +2ω0 light from the circular (0-fold) aperture, respectively.e-f, Same for the structured (3-fold) aperture.These were measured when the simulation domain was centred on z =250λL.
Figure S2: a-b 2D harmonic spectrum obtained for an aperture radius of 1.5 µm and corrugation depths of 0.5 µm and 1.5 µm, respectively.c-d Corresponding harmonic spectra for an aperture radius of 2.5 µm.e-h The 2D spatial profile of the combined 2ω0 generated light for the same four preceding cases.